/* EScheme.h
 *
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 *
 */
#include <src/miscellaneous/OS_basic.h>
#include <src/math/norms.h>
#include <src/math/BoundaryCondition/BoundaryCondition.h>
#include <src/algorithm/Evolution_Algorithm/TimeStepAlgorithm.h>

#if !defined(_EScheme_)
#define _EScheme_
#include <src/algorithm/Evolution_Algorithm/Crank_Nicolson.h>

/*! \brief Interface class for the Entropic scheme Algorithm EScheme.

	This class lists all the variables and virtual functions needed when using the entropic scheme algorithm.
	See the class EScheme for more informations.
*/
template <class TimeTypeOperator>
struct EScheme_Interface
{
	typedef EScheme_Interface													type;
	typedef typename TimeTypeOperator::argument_base_type						argument_base_type;
//	typedef typename TimeTypeOperator::argument_type							u_type;
	typedef Mod_Vector<argument_base_type>										u_type;
	typedef typename Dual<argument_base_type>::type								dual_type;
	OS_STATIC_CHECK(dual_type::Dim>=1);//grr!! always this strange compiler problem.
//	typedef typename TimeTypeOperator::argument_type							vstar_type;
	typedef Mod_Vector<dual_type>												vstar_type;

	EScheme_Interface() : cfl_adjust_(1.), n_(0) {};
//! One call at the beginning of the time step algorithm.
	virtual void init(typename TimeTypeOperator::indice_type n) {};
//! compute vstar_[n]. See class EScheme for more informations
	virtual void vstar(typename TimeTypeOperator::indice_type n) = 0;
//! compute gstar_[n]. See class EScheme for more informations
	virtual void gstar(typename TimeTypeOperator::indice_type n) = 0;
//! return \f$\tau^n\f$. See class EScheme for more informations
	virtual OS_double CFL(typename TimeTypeOperator::indice_type n)  = 0;
//! end criteria for the fixed point algorithm. See class EScheme for more informations
	virtual bool criteria(typename TimeTypeOperator::indice_type n)  = 0;
//! Emergency stop for the fixed point algorithm. See class EScheme for more informations
	virtual bool trigger_adjust_cfl(typename TimeTypeOperator::indice_type n) = 0;
//! This variable stores the solution \f$\{u^n\}_{n\ge 0}\f$.
	u_type u_;
//! This variable stores the solution \f$\{v^{*,k,n+1/2}\}_{n\ge 0}\f$. Convention is vstar_[n] = \f$ v^{*,k,n+1/2}\f$.
	vstar_type vstar_;
//! This variable stores the solution \f$\{g^{*,k,n+1/2}\}_{n\ge 0}\f$. Convention is gstar_[n] = \f$ g^{*,k,n+1/2}\f$.
	vstar_type gstar_;
//! This variable stores the time index n
	typename TimeTypeOperator::indice_type n_;
//! This variable is used for the adaptive time technique.
	typename TimeTypeOperator::base_value_type cfl_adjust_;
//! This variable stores the intermediate time.
	typename TimeTypeOperator::time_type t_12;
//! This variable stores \f$\tau^n\f$.
	typename TimeTypeOperator::time_type delta_t_;
//! Intermediate values \f$u^{*,k,n+1}\f$, \f$\mathcal{O}(u^{*,k,n+1})\f$
	typename TimeTypeOperator::argument_base_type S_k1, Op_S_k1;
//! This variable stores the number of iterations needed for the fixed point algorithm to converge.
	std::size_t nb_iterations_;
};

/*! Entropic schemes of arbitrary orders.
Purposes of these schemes are to compute numerically systems of equations of kind
\f[
    \partial_t u + f (u) = 0,
\f]
where \f$u(t,x) \in \RR^D\f$, \f$x\in \Omega \in \RR^N\f$, and \f$f(u) \in \RR^D \f$ is a (non linear) differential operator.

The system must be writable into the entropy variable form
\f[
    \partial_t u + g(v)= 0
\f]
where \f$v\f$, the entropy variable, is equal to \f$\nabla_u U(u) = v(u)\f$. Here, \f$U(u) \in \RR\f$ is the entropy functional. We assumed thus that
there exists a pair of Entropy / Entropy fluxes \f$(U(u) \in \RR,G(u)\in \RR^N)\f$ such that
\f[
    \partial_t U(u) + \nabla \cdot G(u) \le c(t,x)
\f]
in a weak sense and \f$0\le c\in L^1(\Omega)\f$.

To modelize numerically such systems, let us introduce \f$u^n \sim u(t^n)\f$ the discrete solution at computed time \f$\{t^n\}_{n \ge 0}\f$,
and denote \f$\{\tau^n = t^{n+1}-t^n\}_{n \ge 0}\f$.

The recipee to build a scheme is the following : first choose a \f$(q+1)\f$-time level interpolation \f$u^*(u^q, ..,u^0)\f$ (the maximum order of precision is q+2).
Suppose that :
- This interpolation is consistent with the identity \f$(u^*(u, ..,u) = u)\f$
- The map \f$u^q \rightarrow u^*(u^q, ..,u^0)\f$ is smoothly invertible.
Denote \f$u^{*,n} = u^*(u^n, ..,u^{n-q})\f$.

Now choose the Entropy variable modelling \f$U^*(u^q, ..,u^0)\f$ (\f$U(u^*)\f$ is a classical choice), and denote \f$U^{*,n} = U^*(u^n, ..,u^{n-q})\f$.
It must be
- consistent with the original entropy \f$U(u)\f$ (i.e. \f$U^*(u, ..,u) = U(u)\f$),
- It must define the \f$(q+2)\f$-time entropy variable \f$v^{*,n+1/2}(u^{q+1}, ..,u^0)\f$, that is a variable verifying
\f[
    \delta_t U^{*,n} = \frac{U^{*,n+1}-U^{*,n}}{\tau^n} = v^{*,n+1/2} \cdot \delta_t u^{*,n}
\f]
consistent with the entropy variable : \f$v^{*,n+1/2}(u, ..,u) = v(u)\f$.

The system is then approximated by the numerical scheme
\f[ \label{e-scheme}
    \delta_t u^{*,n} := \frac{u^{*,n+1}-u^{*,n}}{\tau ^n} = - g(v^{*,n+1/2})
\f]
Usually theses schemes are entropy stable under a reasonable hypothesis
: more precisely let \f$E^{*,n} = \sum_i U(u_i^n)\f$. We compute straightforwardly
\f[
    \delta_t E^{*,n} = v^{*,n+1/2} \cdot \delta_t u^{*,n} = v^{*,n+1/2} \cdot g(v^{*,n+1/2})
\f]
Suppose that \f$v^{*,n+1/2} \cdot g(v^{*,n+1/2}) := G(v^{*,n+1/2})\f$ is bounded (i.e.
\f$\sum_i G(v_i^{*,n+1/2}) \le C\f$), then we deduce
\f[
    \delta_t E^{*,n} \le C
\f]
For instance, consider the following special case of an \f$L^2\f$ entropic scheme :
\f[
    u^{*,n} := u^n; U^{*,n} = U(u^n):=\frac{1}{2} |u^n|^2; v^{*,n+1/2} := \frac{1}{2}(u^{n+1}+u^n).
\f]
This setting provides the following numerical scheme :
\f[
    \delta_t u^n = f(v^{*,n+1/2}) = f(u^{n+1/2})
\f]
i.e. Crank Nicolson schemes are a particular case of entropic schemes.


To compute the solution of (\ref{e-scheme}), we use a fixed point algorithm in the Crank Nicolson spirit : define the first iterate
\f[
    u^{0,n+1}:=u^n
\f]
Then, for \f$k>=1\f$, define the entropy variable
\f[
v^{*,k,n+1/2} := v^*(u^{k-1,n+1},u^n,...,u^{n-q})
\f]
compute
\f[
    u^{*,k,n+1} : = u^{*,n} - \tau^n g(v^{*,k,n+1/2});
\f]
deduce the next iterate \f$u^{k,n+1}\f$ from \f$u^{*,k,n+1}\f$.

Loop until a convergence holds.
*/
template <	class TimeTypeOperator = Null >
class EScheme : 
	public TimeStepAlgorithm< TimeTypeOperator >
	, public EScheme_Interface< TimeTypeOperator >
{
public:
	typedef typename Crank_Nicolson<TimeTypeOperator>				Base;
	typedef typename EScheme										algorithm_type;
	virtual ~EScheme(){};
	EScheme() {};

	EScheme(const time_operator_type & Op	) : Base(Op)
	{};
	result_type operator()(const time_type & from, time_type & to, argument_type S_k) 
	{
		nb_iterations_ = 0.;
		delta_t_	= to-from ;
		u_[n_]			= S_k;		// updating computations.
		u_[n_+1]		= S_k;		// updating computations. First iterate equal to identity.
		init(n_);
		vstar(n_);
		gstar(n_);
		delta_t_ =  std::min(delta_t_, cfl_adjust_*CFL(n_));
		OS_DYNAMIC_CHECK(delta_t_ > THRESHOLD,"EScheme : the time step is too low. Check the CFL function.");
		u_[n_+1]    = u_[n_]+delta_t_*gstar_[n_];
			
		t_12		= from + delta_t_ /2. ;

		bool Convergent = false;
		do {
start:				
				nb_iterations_++;

				Op_S_k1 = gstar_[n_];
				S_k1    = u_[n_+1];
				vstar(n_);
				gstar(n_);
				u_[n_+1] = u_[n_] + delta_t_*gstar_[n_];

				if (trigger_adjust_cfl(n_) ) {
					adjust_cfl(from,S_k);
				}
				if ( criteria(n_) ) {
					break;
				}
		} while (true);
		if (nb_iterations_ < 3) {
			if (cfl_adjust_ < 1) //we should never overpass the cfl condition, that ensures stability !
				cfl_adjust_ *= 2.; // multiply CFL constant
		}
		to = from + delta_t_;
		++n_;
		return u_[n_];
	}
// a standard implementation
	virtual bool trigger_adjust_cfl(indice_type n)  {
			return (nb_iterations_ > 10 );
	};

	virtual void update(){
		typedef u_type::indice_list_type indice_list_type;
		if (u_[n_].capacity() != get_state_variable()->capacity()) {
			indice_list_type indx = u_.size_list();
			indice_list_type indx_source = get_state_variable()->size_list();
			std::copy(indx_source.begin(),indx_source.end(),indx.begin());
			u_.resize(indx);
			gstar_.resize(indx);
			vstar_.resize(indx);
			u_[n_] = *get_state_variable();
		};
		Op_.set_state_variable(get_state_variable());
	}; 
	static std::string get_name(){return " E Scheme delta_t u = ";};
protected :
	virtual void adjust_cfl(const time_type & from, argument_type S_k){
		OS_DYNAMIC_CHECK(cfl_adjust_ > THRESHOLD,"EScheme : no convergence of the fixed point algorithm. Check the convergence criteria..");
		cfl_adjust_ /= 2.; // divide CFL constant
		delta_t_ /= 2.; 
		t_12		= from + delta_t_ / 2.;
		u_[n_+1] = u_[n_];
		vstar(n_);
		gstar(n_);
		u_[n_+1] = u_[n_]+delta_t_* gstar_[n_];
		nb_iterations_ = 0;
	}
};

#endif
